# DGP for Signorino Selection (private information)
rm(list = ls())

library(MSBVAR)
library(foreign)  					# Packages for data
library(mnormt)						# Packages for MVNormal
library(mvtnorm)					# Packages for MVNormal

#################################################
# The Code: Structural Estimator NO Selection
#################################################

struct <- function(x11,x14,x24,y,beta){
	x11=cbind(1,as.matrix(x11))
	x14=as.matrix(x14)
	x24=as.matrix(x24)
    x24 <- cbind(1,as.matrix(x24))
	n <- dim(x11)[1]
	n11 <- dim(x11)[2]
	n14 <- dim(x14)[2]
	n24 <- dim(x24)[2]
	b11 <- as.matrix(beta[1:n11])
	b13 <- as.matrix(beta[(n11+1)])
	b14 <- as.matrix(beta[(n11+2):(n11+2+n14-1)])
	b24 <- as.matrix(beta[-c(1:(n11+2+n14-1))])
	u24 <- x24%*%b24
	#p2 <- pnorm((u24)/sqrt(2))
	p2 <- pnorm((u24))
	u11 <- x11%*%b11
	u13 <- b13*rep(1,n)
	u14 <- x14%*%b14
	u12 <- p2*u14 + (1-p2)*u13
	#p1 <- pnorm((p2*u14+(1-p2)*u13-u11)/sqrt(p2^2+(1-p2)^2+1))
	p1 <- pnorm((p2*u14+(1-p2)*u13-u11))
	vec1 <- which(y[,1]==0)
	vec3 <- which(y[,2]==0)
	vec4 <- which(y[,2]==1)
	chunk1 <- log(1-p1)
	chunk3 <- log(p1*(1-p2))
	chunk4 <- log(p1*p2)
	ll <- sum(chunk1[vec1]) + sum(chunk3[vec3]) + sum(chunk4[vec4]) 
	return(-ll)	
	}





#################################################
# The Code: Selection Estimator NO Structure WITH NORMAL SUBSCRIPTS
#################################################

# notation reasonable

selec12 <- function(x1,x2,y,beta){
	# NOTE: parametrization DOES NOT follow Greene/Wawro!!!!
	x1 <- cbind(1,x1)										# adding a constant 
	x2 <- cbind(1,x2)										# adding a constant
	n <- dim(x1)[1]; a <- dim(x1)[2]; b <- dim(x2)[2]		# prepare to slice the coef vector
	C1 <- a+b; a1 <- a+1; raus <- c(1:C1)					# prepare to slice the coef vector
	para <- beta[-raus]										# Reparametrizing Rho
	rho <- 2*(1/(1+exp(-para)))-1
	b1 <- beta[1:a]											# The Beta's
	b2 <- beta[a1:C1]
	xb1 <- x1%*%b1; xb2 <- x2%*%b2							# xb's
	#part1 <- cbind(-xb1,xb2); part2 <- cbind(xb1,xb2)		# combining for mvnorm
	part1 <- cbind(-xb2,xb1); part2 <- cbind(xb2,xb1)		# combining for mvnorm
	covar1 <- matrix(c(1,-rho,-rho,1),2,2)
	covar2 <- matrix(c(1,rho,rho,1),2,2)
	prob1 <- apply(part1,1,pmnorm, mean=rep(0,2),varcov=covar1)	# y_1=1; y_2=0
	prob2 <- apply(part2,1,pmnorm, mean=rep(0,2),varcov=covar2)	# y_1=1; y_2=1
	prob3 <- 1-pnorm(xb1)										# y_1=0; y_2=?
	Y1 <- y[,1]; Y2 <- y[,2]
	vec3 <- which(Y1==0); vec2 <- which(Y2==1); vec1 <- which(Y2==0) # which obs belong where?
	chunk1 <- sum(log(prob1[vec1]))								# for each y group, the 
	chunk2 <- sum(log(prob2[vec2]))								# partial LL contrinution
	chunk3 <- sum(log(prob3[vec3]))								# it makes
	ll <- chunk1 + chunk2 + chunk3 
	return(-ll)
	}





#################################################
# The Code: Selection A N D  Structure Estimator
#################################################


selstruc <- function(x11,x14,x24,y,beta){
	# NOTE: parametrization is now intuitive (first 1, then 2) unlike Greene/Wawro
	x11=cbind(1,as.matrix(x11))
	x14=as.matrix(x14)
	x24=as.matrix(x24)
    x24 <- cbind(1,as.matrix(x24))
	n <- dim(x11)[1]
	n11 <- dim(x11)[2]
	n14 <- dim(x14)[2]
	n24 <- dim(x24)[2]
	b11 <- as.matrix(beta[1:n11])
	b13 <- as.matrix(beta[(n11+1)])
	b14 <- as.matrix(beta[(n11+2):(n11+2+n14-1)])
	b24 <- as.matrix(beta[(n11+2+n14):(n11+2+n14+n24-1)])
	para <- as.matrix(beta[-c(1:(n11+2+n14+n24-1))])
	rho <- 2*(1/(1+exp(-para)))-1
	xb2.sp <- x24%*%b24												# latent D.utility for player 2 for a_4
	p.p <- pnorm(xb2.sp/(sqrt(1)))									# p(Y4=1|Y1=0)
	u1.34 <- -x11%*%b11 + p.p*(x14%*%b14) + (1-p.p)*rep(b13,n)		# latent D.utility for player 1 for a_2
	part1 <- cbind(-xb2.sp,u1.34); part2 <- cbind(xb2.sp,u1.34)		# combining for mvnorm
	covar1 <- matrix(c(1,-rho,-rho,1),2,2)				# negative rho
	covar2 <- matrix(c(1,rho,rho,1),2,2)				# positive rho
	prob1 <- apply(part1,1,pmnorm, mean=rep(0,2),varcov=covar1)	# y_1=1; y_2=0
	prob2 <- apply(part2,1,pmnorm, mean=rep(0,2),varcov=covar2)	# y_1=1; y_2=1
	prob3 <- 1-pnorm(u1.34)										# y_1=0; y_2=?
	#Y1 <- y[,1]; Y2 <- y[,2]
	Y <- y
	#vec3 <- which(Y1==0); vec2 <- which(Y2==1); vec1 <- which(Y2==0) # which obs belong where?
	vec3 <- which(Y==1); vec2 <- which(Y==4); vec1 <- which(Y==3) # which obs belong where?
	chunk1 <- sum(log(prob1[vec1]))								# for each y group, the 
	chunk2 <- sum(log(prob2[vec2]))								# partial LL contrinution
	chunk3 <- sum(log(prob3[vec3]))								# it makes
	ll <- chunk1 + chunk2 + chunk3 
	return(-ll)
	}


#### 	D	A	T	A




start.time <-  Sys.time()
m <- 100
N <- 1000
#set.seed(020005001)
rr <- rep(NA,m)
parameters.selec <- matrix(NA,m,7)
parameters.struc <- matrix(NA,m,7)
parameters.selStruc <- matrix(NA,m,8)
ses.selStruc <- matrix(NA,m,8)
parameters.prob.1 <- matrix(NA,m,3)
parameters.prob.2 <- matrix(NA,m,3)
RHO <- 0
VAR.x <- 1
# conditions: 
#rho -0.9, -0.5, 0.0, 0.5, 0.9 
# VAR.x= 0.5; 1.0; 2.0; 10
# N = 200 / 500 / 1000

# var-ratio=(\sigma^2_{eps1}+\sigma^2_{eps2})/(2\sigma^2_x)
# conditions: var-ratio = 0.1, 0.5, 1, 2

for (j in 1:m){
varcov <- matrix(c(1,RHO, RHO, 1), 2,2)
e <- rmvnorm(N, c(0,0), varcov)
X1 <- cbind(rnorm(N,0,VAR.x), rnorm(N,0,VAR.x))
X2 <- cbind(rnorm(N,0,VAR.x), X1[,2])
e1 <- e[,1]
e2 <- e[,2]

# create second node actions
Y2 <- 0+ X2[,1] + X2[,2] + e2
for (i in 1:N){
	ifelse(Y2[i]<0, Y2[i] <- 0, Y2[i] <- 1)	
	}
Y2.star <-  X2[,1] + X2[,2]
cond.mean <- 0 + 1/1*(RHO)*(e1-0) # this deviates from wikip. b/c 1/rho
cond.var <- (1-RHO^2)
e2.hat <- rnorm(N,mean=cond.mean,sd=sqrt(cond.var))
p.expec <- pnorm(Y2.star + e2.hat)

# create first node actions
Y1 <- 0 - X1[,1] + p.expec * X1[,2] + e1
for (i in 1:N){
	ifelse(Y1[i]<0, Y1[i] <- 0, Y1[i] <- 1)	
	}
# eliminate unobserved outcomes
for (i in 1:N){
	ifelse(Y1[i]==0, Y2[i] <- NA, Y2[i] <- Y2[i])
	}
# create outcome matrix
Y <- cbind(Y1,Y2)
YY <- rep(NA, dim(Y)[1])
for (i in 1:length(YY)){
	ifelse(Y[i,1]==0, YY[i] <- 1, ifelse(Y[i,2]==1,YY[i]<-4,YY[i]<-3 ))
}


# Actual correlation and simple GLM's
rr[j] <- cor(e)[2,1]
out.bias.1 <- glm(Y1 ~ X1, family=binomial(link="probit"))
out.bias.2 <- glm(Y2 ~ X2, family=binomial(link="probit"))
parameters.prob.1[j,] <- summary(out.bias.1)$coef[,1]
parameters.prob.2[j,] <- summary(out.bias.2)$coef[,1]



# Estimation Structural
#startval <- c(summary(out.bias.1)$coef[,1],summary(out.bias.1)$coef[,2])
startval <- rep(-0.1,7)
out.struc <- optim(startval, struct, x11=X1[,1],x14=X1[,2], x24=X2, y=Y, hessian=TRUE, 	method="BFGS", control=list(trace=5, REPORT=4))
parameters.struc[j,] <- out.struc$par


# Estimation Select
startval <- c(summary(out.bias.1)$coef[,1],summary(out.bias.1)$coef[,2],0)
out.selec <- optim(startval, selec12, x1=X1, x2=X2,y=Y, hessian=TRUE, 	method="BFGS", control=list(trace=5, REPORT=4))
parameters.selec[j,] <- out.selec$par
parameters.selec[j,7] <- 2*(1/(1+exp(-out.selec$par[7])))-1
#ses.selec[j,] <- sqrt(diag(try(solve(out.selec$hessian))))

# Estimation Select&Struct
#startval <- c(summary(out.bias.1)$coef[,1],summary(out.bias.1)$coef[,2],-.1)
startval <- rep(-0.1,8)
out.selStruc <- optim(startval, selstruc, x11=X1[,1], x14=X1[,2], x24=X2,y=YY, hessian=TRUE, 	method="BFGS", control=list(trace=5, REPORT=4))
parameters.selStruc[j,] <- out.selStruc$par
#parameters.selStruc[j,8] <- 2*(1/(1+exp(-out.selStruc$par[8])))-1
ses.selStruc[j,] <- try(sqrt(diag(solve(out.selStruc$hessian))), silent=TRUE)



 if (j %% 1 == 0) print(paste("Finished MC Simulation iteration",j))
}

end.time <-  Sys.time()
cat(paste("Job finished at:",end.time))
dur <- end.time-start.time
cat(paste("Job took:",dur))


stt <- format(start.time, "%X")
fit <- format(end.time, "%X")

parameters.struc
parameters.selec
parameters.selStruc

colMeans(parameters.struc)
colMeans(parameters.selec)
colMeans(parameters.selStruc)

bias.struc <- (c(0,1,0,1,0,1,1) - colMeans(parameters.struc))^2
bias.selStruc <- (c(0,1,0,1,0,1,1) - colMeans(parameters.selStruc)[-8])^2
bias.selec <- (c(0,1,1,0,1,1) - colMeans(parameters.selec)[-7])^2
bias.selec <- c(bias.selec[1:2],10,bias.selec[3:6])


timeR <- c(stt,fit)

stuff <- list(param.struc=parameters.struc, param.selec=parameters.selec, param.selStruc=parameters.selStruc,SE.selStruc=ses.selStruc, corr.rho=rr, RunTime=timeR, bias.struc=bias.struc, bias.selStruc=bias.selStruc, bias.selec=bias.selec)

save(stuff,file= "/Users/lucas/Dropbox/Methods Paper/Selection/Results revised MC/ResultSimulation_00.RData")